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Abstract: We consider optimal control problems of elliptic PDEs on hypersurfaces T in K n for 
n = 2, 3. The leading part of the PDE is given by the Laplace-Beltrami operator, which is discretized by 
finite elements on a polyhedral approximation of T. The discrete optimal control problem is formulated 
on the approximating surface and is solved numerically with a semi-smooth Newton algorithm. We 
derive optimal a priori error estimates for problems including control constraints and provide numerical 
examples confirming our analytical findings. 
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1 Introduction 

We are interested in the numerical treatment of the following linear-quadratic optimal control 
problem on a re-dimensional, sufficiently smooth hypersurface T C M n+1 , n = 1,2. 

77 \ 1 ll l|2 , a i| ||2 

mm J(u,y) = - \\y — z r2/ri H — ii 

subject to u £ U a d and (1-1) 
J VryVr^ + cy^dT = J u<p dT , V<^ G H 1 (T) 

with U ai i = {v G L 2 (T) \ a < v < 6} , a < b G M . For simplicity we will assume T to be 
compact and c = 1. In section [4] we briefly investigate the case c = 0, in section [5] we give an 
example on a surface with boundary. 



Problem (1.1) may serve as a mathematical model for the optimal distribution of surfactants 
on a biomembrane V with regard to achieving a prescribed desired concentration z of a 
quantity y. 



It follows by standard arguments that (1.1) admits a unique solution u G U a d with unique 
associated state y = y{u) G H 2 (T). 



Our numerical approach uses variational discretization applied to (1.1), see [Hin05j and 



|HPUU 09j . on a discrete surface F h approximating T. The discretization of the state equation 



in (1.1) is achieved by the finite element method proposed in |Dzi88| , where a priori error 
estimates for finite element approximations of the Poisson problem for the Laplace-Beltrami 
operator are provided. Let us mention that uniform estimates are presented in |Dcm09] , and 
steps towards a posteriori error control for elliptic PDEs on surfaces are taken by Demlow and 
Dziuk in [DD07]. For alternative approaches for the discretization of the state equation by 
finite elements see the work of Burger |Bur08| . Finite element methods on moving surfaces 
are developed by Dziuk and Elliott in [DE07J. To the best of the authors knowledge, the 
present paper contains the first attempt to treat optimal control problems on surfaces. 

We assume that T is of class C 2 with unit normal field v. As an embedded, compact hy- 
persurface in M n+1 it is orientable and hence the zero level set of a signed distance function 
\d(x)\ = dist(x,r). We assume w.l.o.g. Vd{x) = v(x) for x G T. Further, there exists an 
neighborhood M C M n+1 of T, such that d is also of class C 2 on J\f and the projection 
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a:Af-tT, a(x) = x - d(x)X7d(x) (1.2) 
is unique, see e.g. |GT981 Lemma 14.16]. Note that Vd(x) = u{a{xj). 

Using a we can extend any function (j) : T — > M to N as <fi(x) = <p{a{x)). This allows us to 
represent the surface gradient in global exterior coordinates Vr0 = {I — vv T )V(f), with the 
euclidean projection (/ — vv T ) onto the tangential space of V. 

We use the Laplace-Beltrami operator Ap = Vr • Vr in its weak form i.e. Ar : H l {T) — > 

y ^ - ^V r yV r (-)dr G H\r)* . 

Let S denote the prolongated restricted solution operator of the state equation 

S : L 2 {T) -»• L 2 (r) , u^y - A r y + cy = u , 

which is compact and constitutes a linear homeomorphism onto H 2 (T), see [Dzi88[ 1. Theo- 
rem]. 

By standard arguments we get the following necessary (and here also sufficient) conditions 
for optimality of u G U a d 

(V u J(u,y(u)),v - u) L 2(r) = {au + S*(Su - z),v - u) L 2 (r) > \/v£U ad , (1.3) 



We rewrite (1.3) as 

u = P Uad (-~S*(Su-z) \ , (1.4) 
where Pu ad denotes the L 2 -orthogonal projection onto U a d- 

2 Discretization 

We now discretize (1.1) using an approximation T h to T which is globally of class C 0,1 . 



Following Dziuk, we consider polyhedral T h = Uie/ h ^-h consisting of triangles with corners 
on r, whose maximum diameter is denoted by h. With FEM error bounds in mind we assume 
the family of triangulations T h to be regular in the usual sense that the angles of all triangles 
are bounded away from zero uniformly in h. 

We assume for T h that a(F h ) = T, with a from (1.2). For small h > the projection a also 



is injective on Y . In order to compare functions defined on T h with functions on Y we use a 
to lift a function y G L 2 (T h ) to T 

y\a{x)) = y(x) Vx G T h , 

and for y G L 2 (T) and sufficiently small h > we define the inverse lift 

Vi{x) = y{a{x)) \/x£T h . 

For small mesh parameters h the lift operation (•)/ : L 2 (T) — > L 2 (T h ) defines a linear homeo- 
morphism with inverse (•)'. Moreover, there exists c- m t > such that 

1 - Cmt/i 2 < ll(-)zll£(L 2 (r),L 2 (r' i ))> ll(')'ll£(L 2 (r' i ),L 2 (r)) — 1 + c mth 2 , (2.1) 
as the following lemma shows. 
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Lemma and Definition 2.1. Denote by the Jacobian of a\^h : T h — > V, i.e. = 
|det(M)| where M G M. nxn represents the Derivative da(x) : T x T h — > T a ^T with respect to 
arbitrary orthonormal bases of the respective tangential space. For small h > there holds 



sup 
r 



1 



dr 



dT h 



< c int h , 



Now let ^rp- denote |det(M 1 )| ; so that by the change of variable formula 



[ vi dT h - [ v dr = / 



< Cint/l 2 ||v| 



Proof, see [DE07, Lemma 5.1] 



Problem (1.1) is approximated by the following sequence of optimal control problems 

1„ 



min J( u , y) 



a . 



IV ~ z i\\L 2 (r h ) + 2 \\ u \\~L 2 (r h ) 



subject to u G U^ d and 
V = S h u , 



(2.2) 



with U^ d = {i)E! L 2 (T h ) \ a < v < 6}, i.e. the mesh parameter h enters into U a d only through 
T h . Problem ( |2.2[ ) may be regarded as the extension of variational discretization introduced 
in [Hin05| to optimal control problems on surfaces. 

In |Dzi88j it is explained, how to implement a discrete solution operator Sh '■ L 2 (T h ) — > 
L 2 (T h ), such that 

11(0 Sh(-)i - S\\ C (L2(r),L 2 (r)) < Cfe^ 2 , (2.3) 

which we will use throughout this paper. See in partikular [Dzi88, Equation (6)] and [Dzi88, 
7. Lemma]. For the convenience of the reader we briefly sketch the method. Consider the 
space 



V h 



<pec°(r h 



ViG I h : V \ Tl ^V l {r h )]dH l (T h ) 



of piecewise linear, globally continuous functions on V . For some u G L 2 (r), to compute 
Vh = (-y s h{-)iu solve 

/ Vrh^Vph^i + cy h <fi dT h = mpi dT h , V<p€V h 
Jr h Jr h 

for yh G Vh- We choose L 2 {T h ) as control space, because in general we cannot evaluate f r v dT 
exactly, whereas the expression f rh vi dT h for piecewise polynomials vi can be computed up 
to machine accuracy. Also, the operator Sh is self-adjoint, while ((O'S/iCOO* = (-)i* Sh(-) 1 * is 
not. The adjoint operators of (•); and (•) have the shapes 



to G L 2 (T h ) : ((■),)*« 



dT h 
~dY 



\/v G L 2 (T) : ((•)')*« 



dr 

df^ 



v, 



(2.4) 



hence evaluating (-)z* and (•)'* requires knowledge of the Jacobians ^jp- and which may 
not be known analytically. 



3 



Similar to (1.1), problem (2.2) possesses a unique solution Uh G Uff d which satisfies 



a 



~Ph(Uh) 



(2.5) 



Here P V h : L 2 (T h ) — »■ UK is the L 2 (r ft )-orthogonal projection onto UK and for v G L 2 (r fe ) 
the adjoint state is Ph(v) = Sl(ShV — zi) G -ff 1 (r /l ). 

Observe that the projections Pu ad and Pjjh coincide with the point-wise projection P[ a ,b] on 
T and r , respectively, and hence 



(2.6) 



for any v G L 2 (r). 

Let us now investigate the relation between the optimal control problems (1.1) and (2.2). 



Theorem 2.2 (Order of Convergence). Let u G L (T), Uh G L (T ) be the solutions of (1.1) 
and (2.2) ; respectively. Then for sufficiently small h > there holds 

1 + Cint/l 2 fl 



n|| "/< ~ a ilz, 2 (r) + \\y"h - y\ 



< 



^(n-i-c^u 



2 

L3(T) 



+ 



5 U 



L2(r) 



(2.7) 



mi/i y = 5u and yh = S^Uh- 

Proof. From (2.6) it follows that the projection of — (aP(u)), onto U% d is ui 

1 



1i/ = Prjh 



a 



-P( u )l 



which we insert into the necessary condition of (2.2). This gives 

(au h +p h (u h ),ui - u h ) L 2( V h) > 0. 

On the other hand U\ is the L 2 (r /l )-orthogonal projection of —^p(u)i, thus 

1 

( p(u)i -u h u h - ui) L 2 (V h) < . 

Adding these inequalities yields 

a\\ui - Uh\\ L 2( T h) <({Ph{uh) ~p(u)i) ,ui - u h ) L 2(ph) 

={Ph(uh) - S* h {y - z)i,ui - u h ) L 2(Y h) + {S^iy - z)i - p(u)i,ui - u h ) L 2( Th) 

The first addend is estimated via 

{Ph(uh) - Sl(y - z)i,ui - u h ) L 2 {V h) = (y h - yi,S h ui - yh)i 2 (r h ) 

2 



\yh - yi\\ L 2(rh) + (yh ~ Vh Shn - yi) L 2 (r h ) 



< 



1 



1 



yh - yi\\ L 2(rh) + »\\ s hUi - yi\\ L 2 {rh) 
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The second addend satisfies 

{Sh(y - z)i-p(u)i,ui - u h ) L 2 {vh) < - u h \\ 2 L2{rh) + ^\\s*h(y - z)i - p{ u )i\\ 2 L 2 {vh) ■ 

Together this yields 

"IN - u h \\ 2 L2{vh) + \\y h - yz||| 2 ( r fc) < ~W S h(y ~ z )i ~ P( u )i\\ 2 Li(Th) + \\ShUi - yi\\ 2 L 2 {rh) 



The claim follows using (2.1) for sufficiently small h > 



Because both S and Sh are self-adjoint, quadratic convergence follows directly from (2.7). 
For operators that are not self-adjoint one can use 



■)i*Sh('Y - S'*ll,c(L 2 (r),L2(r)) < C FE h z . 



(2.8 



which is a consequence of (2.3). Equation (2.4) and Lemma |2.l| imply 

- (■) l \\c(L2(r»),L*(r)) < Cint/i 2 , ||((-) Z r - (■)i\\c(LHr),m^)) < c int h 2 . (2.9) 



Combine ( |2.7| ) with (2.8) and (2.9) to proof quadratic convergence for arbitrary linear elliptic 
state equations. 



3 Implementation 



In order to solve ( |2.5| ) numerically, we proceed as in |Hin05| using the finite element techniques 
for PDEs on surfaces developed in [Dzi88] combined with the semi-smooth Newton techniques 
from [HIK03J and [Ulb03] applied to the equation 



G h (u h ) = \ U h - P[ atb ] (-^Ph(Uh) 







(3-1) 



Since the operator ph continuously maps v 6 L 2 (T h ) into if 1 (r /l ), Equation (3.1) is semis- 
mooth and thus is amenable to a semismooth Newton method. The generalized derivative of 
G h is given by 

X 



DG h (u) = lI+^S* h S h ) , 



o 



where x '■ ^ h — > {0,1} denotes the indicator function of the inactive set X(— -ph(u)) 

{ 7 er h |a<-i Pfc («)[ 7 ]<6} 



X 



1 onZ(-WtO) cT h 



elsewhere on T 



I, 



which we use both as a function and as the operator \ : L 2 (T h ) — > L' 2 (T h ) defined as the 



-2/rh\ 



point-wise multiplication with the function \- A step semi-smooth Newton method for (3.1) 
then reads 



I+^S* h S h )u" 
a 



-G h (u) +DG h (u)u = Pr a6 i (--p h (u)) +-S* h S h u. 

1 ' \ a J a 



Given u the next iterate u + is computed by performing three steps 
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1. Set ((1 - x) = ((1 - x)P[a,b] {-^Ph(u) + rn)) [7], which is either a or 6, depend- 
ing on 7 G T h . 

2. Solve 

(/ + ^S* h S h ) X u + = \ [Slzi - S* h S h (1 - X ) u + ) 
for x u+ by CG iteration over L 2 (I(— ^ph(u)). 

3. Set u" 1 " = + (1 — x) u+ ■ 
Details can be found in [H V11] . 



4 The case c = 

In this section we investigate the case c = which corresponds to a stationary, purely diffusion 
driven process. Since T has no boundary, in this case total mass must be conserved, i.e. the 
state equation admits a solution only for controls with mean value zero. For such a control 
the state is uniquely determined up to a constant. Thus the admissible set U a d has to be 
changed to 



U ad 



{v G L 2 (T) \ a<v< b} DLg(r), where Lg(T) := iu G L 2 (T) J vdT = j 



and a < < b. Problem (1.1) then admits a unique solution (k, y) and there holds J" r ydr = 
J* r z dr. W.l.o.g we assume J r z dT = and therefore only need to consider states with mean 
value zero. The state equation now reads y = Su with the solution operator S : Lq(F) — > 
Lq(T) of the equation — Ary = u, J r ydT = 0. 

Using the injection Lq(T) A L 2 (T), S is prolongated as an operator S : L 2 (T) — > L 2 (T) by 
S = 1S1*. The adjoint 1* : L 2 (T) — > Ll(T) of 1 is the L 2 -orthogonal projection onto 



The unique solution of (1.1) is again characterized by (1.4), where the orthogonal projection 
now takes the form 



with m G M chosen such that 



F u ad (v) = P[ a ,b] + m) 
J P [a>b] (v + m) dr = o. 



If for v G L 2 (T) the inactive set T(v + rn) = {7 G T \ a < v[y] + m < b} is non-empty, the 
constant m = m(v) is uniquely determined by v G L 2 (T). Hence, the solution u G U a d satisfies 



u 



[0,6] 



—p(u) + m ( — p(u) 
a \ a 



withp(n) = S*(Su—i*z) G H 2 (T) denoting the adjoint state and m(— — p(u)) G M is implicitly 
given by J r udT = 0. Note that 1*1 is the identity on Lq(T). 

In we now replace U^ d by [/ a h d = {v G L 2 (r h ) | a < v < b} n Lg(r fe ). Similar as in glj), 
the unique solution Uh then satisfies 



a 



-Phiuh) 



a,b] I Ph{uh) +m h 



a 



a 



~Ph(Uh) 



(4.1) 
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with Ph(vh) = S^ShVh — i* h zi) G H l {V h ) and rnh(—^ph(uh)) G K the unique constant such 
that J rh n/j dr' 1 = 0. Note that (—^Ph(uh)) is semi-smooth with respect to Uh and thus 
Equation (4.1) is amenable to a semi-smooth Newton method. 

The discretization error between the problems (2.2) and (1.1) now decomposes into two 
components, one introduced by the discretization of U a d through the discretization of the 
surface, the other by discretization of S. 

For the first error we need to investigate the relation between Pjjh (u) and Pu ad ( u ), which 

ad 

is now slightly more involved than in (2.6). 

Lemma 4.1. Let h > be sufficiently small. There exists a constant C m > depending only 
on T, \a\ and \b\ such that for all v G L 2 (T) with fx(v+m(v)) dr > there holds 

\m h (vi) - m(v)\ < — — h 2 



X(v+m(v)) 



dr 



Proof. For v G L 2 (T), e > choose 5 > and h > so small that the set 
Xl = | 7 G T h | a + 5 < vi(i) + m(v) < b - <$} . 

satisfies f xS dT h (l + e) > $x(v+m{v)) ^ ^ s eas y ^° snow that hence m^vi) is unique. Set 
C = Ci nt max(|a|, \b\) f r dT. Decreasing h further if necessary ensures 



For x G R let 



M*{x) = [ P M] («, + X) dT h . 



Since J* r P[ a ,6] ( v + m(v)) dT = 0, Lemma 2.1 yields 

\M*{m{v))\ < c int ||P M (v + m{v)) \\ Ll{T) h 2 < Ch 2 . 
Let us assume w.l.o.g. —Ch 2 < M^(m(v)) < 0. Then 

I Ch 2 ' 



M£ m(v) + 



implies < m(v) — muiv) < f , rh < j — - 

Jl* ai JX(v 

to x. This proves the claim. 



r ^ j >M^(m(v)) + Ch 2 >0 

( 1+e ^ c — h 2 , since M^{x) is continuous with respect 



(v-\-m(v)) 



dr 



Because 



P r/ fc. (vi 



F U ad 0) = P[ a ,6] (« + - P[ a , 6 ] (v + m(«)) , 



we get the following corollary. 

"orollary 4.2. L 

€L 2 (T) «^/ x(w+to(i0) 

Si« (t;i) ) , ~ p ^ {v) 



Corollary 4.2. Lei h > be sufficiently small and Cm as in Lemma J^.l. For any fixed 
v G L 2 (T) with f-r,.. , dr > we have 



/rdr 



L»(T) 



X(i>+m(i>)) 



dr 
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Note that since for u E L 2 (F) the adjoint p(u) is a continuous function on T, the corollary is 
applicable for v = ——p(u). 



The following theorem can be proofed along the lines of Theorem 2.2 

Theorem 4.3. Let u G L 2 (T), Uh G L 2 (T h ) be the solutions of (1.1) and (2.2), respectively, 



in the case c = 0. Let Uh = [Pjjh {—^p(u)i) ) • Then there holds for e > and < h < h 



«K - «hllia(r) + IW - y|lL»(r) ^ ( x + e ) Q || ((O'-SaCOi - (y - *) 



2 

L2(r) 



+ 



(0 l 5&(-)jfifc - y 



L2(r) 



Using Corollary 4.2 we conclude from the theorem 



\ u h ~ HU 2 (r) <C 



)(y-z) 


i 


+ -F= 




L 2 (r) v« 



L2(r) 



+ 



A | \\S\\c{L^T)^(T)) 



c m J/ r dr/i 2 



Jz(-ip(«)+m(-ip(«))) 



dr 



the latter part of which is the error introduced by the discretization of U a d- Hence one has 
/t 2 -convergence of the optimal controls. 

5 Numerical Examples 

The figures show some selected Newton steps u + . Note that jumps of the color-coded function 
values are well observable along the border between active and inactive set. For all examples 
Newton's method is initialized with uq = 0. 

The meshes are generated from a macro triangulation through congruent refinement, new 
nodes are projected onto the surface T. The maximal edge length h in the triangulation is 
not exactly halved in each refinement, but up to an error of order 0(h 2 ). Therefore we just 
compute our estimated order of convergence (EOC) according to 

™^ ln W u hi-i ~ u i\\ L 2(r h i-i) - ln ll n ^ - u i\\^(r h i) 
EOa = ~ H2T ' 

For different refinement levels, the tables show L 2 -errors, the corresponding EOC and the 
number of Newton iterations before the desired accuracy of 10~ 6 is reached. 
It was shown in |HU04| . under certain assumptions on the behaviour of — ^-p(u), that the 
undamped Newton Iteration is mesh-intdependent. These assumptions are met by all our 
examples, since the surface gradient of ——p(u) is bounded away from zero along the border 
of the inactive set. Moreover, the displayed number of Newton-Iterations suggests mesh- 
independence of the semi-smooth Newton method. 

Example 5.1 (Sphere I). We consider the problem 



min J(u, y) subiect to — Arw + y 



u — r, 



-1< u < 1 



(5.1) 
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Step 3 



Step 6 




Figure 1: Selected full Steps u + computed for Example 5.1 on the twice refined sphere 



reg. refs. 





1 2 


3 4 


5 


L2-error 


5.8925e-01 


1.4299e-01 3.5120e-02 


8.7123e-03 2.2057e-03 


5.4855e-04 


EOC 




2.0430 2.0255 


2.0112 1.9818 


2.0075 


# Steps 


6 


6 6 


6 6 


6 




Table 1: L 


2 -error, EOC and number 


of iterations for Example 


5.1 


with r the unit sphere 


in R 3 and a = 1.5 • 10~ 6 . 


We choose z = 52axs(xf 


- x|) , to obt 



the solution 

u = r = min (l, max ( — 1, 4:X^{x\ — a?!))) 



of (5.1) 



Example 5.2. Let T = {{x\, x 2 , x^) T G M 3 | x 3 = x±x 2 A x±,x 2 G (0, 1) } and a = 10 3 . For 
min J(u, y) subject to — Ary = u — r, y = on dT — 0.5 < u < 0.5 

we get 

u = r = max ( — 0.5, min (0.5, sin(-7rx) sin(7ry))) 
by proper choice of z (via symbolic differentiation). 



Example 5.2 although c = 0, is also covered by the theory in Sections T]|3 as by the Dirichlet 
boundary conditions the state equation remains uniquely solvable for u E L 2 (T). In the last 
two examples we apply the variational discretization to optimization problems, that involve 
zero-mean- value constraints as in Section HI 
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reg. refs. 





1 


2 


3 


4 


5 


L2-error 


3.5319e-01 


6.6120e-02 


1.5904e-02 


3.6357e-03 


8.8597e-04 


2.1769e-04 


EOC 




2.4173 


2.0557 


2.1291 


2.0369 


2.0250 


# Steps 


11 


12 


12 


11 


13 


12 



Table 2: L -error, EOC and number of iterations for Example 5.2 



Step 1 



Step 3 



Step 4 



StepS 




Figure 3: Selected full Steps u + computed for Example 5.3 on once refined sphere 



reg. refs. 


1 2 


3 4 


5 


L2-error 


6.7223e-01 1.6646e-01 4.3348e-02 


1.1083e-02 2.7879e-03 


6.9832e-04 


EOC 


2.0138 1.9412 


1.9677 1.9911 


1.9972 


# Steps 


8 8 7 


7 6 


6 




Table 3: L 2 -error, EOC and number 


of iterations for Example 


5.3 



Step 2 



Step 3 



Step 6 




Figure 4: Selected full Steps u + computed for Example 5.4 on the once refined torus. 



reg. refs. 





1 


2 


3 


4 


5 


L2-error 


3.4603e-01 


9.8016e-02 


2.6178e-02 


6.6283e-03 


1.6680e-03 


4.1889e-04 


EOC 




1.8198e+00 


1.9047e+00 


1.9816e+00 


1.9905e+00 


1.9935e+00 


# Steps 


9 


3 


3 


3 


2 


2 



Table 4: L -error, EOC and number of iterations for Example 5.4 
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Example 5.3 (Sphere II). We consider 

min J(u, v) subject to — Ary = u , — 1 < u < 1 , / -yaT= / uaT = 0, 

with r the unit sphere in M 3 . Set a = 10 -3 and 



z(xi,x 2 ,x 3 ) = Aax 3 + < 



ln(x 3 + 1) + C , if 0.5 < x 3 

x 3 — |arctanh(x3) , if —0.5 < x 3 < 0.5 , 
-C-ln(l-x 3 ), if x 3 < -0.5 



where C is chosen for z to be continuous. The solution according to these parameters is 

u = min (l, max ( — 1, 2x3)) . 
Example 5.4 (Torus). Let a = 10~ 3 and 



r = < 



(xi,X2,x 3 ) T G 



^3+ V X 1 +X 2- 1 



2 1 



2 

the 2- Torus embedded in M 3 . By symbolic differentiation we compute z, such that 

min J(u, y) subject to — A-py = u — r, — 1 < u < 1 , ydT = udT 
?),yem(r) Jr Jr 



nSL 2 (r) 

is solved by 



it = r 



max ( — 1, niin (l, hxyz)} 



As the presented tables clearly demonstrate, the examples show the expected convergence 
behaviour. 
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